#include <vector>
#include <iostream>
#include <math.h>
#include <iomanip>
using namespace std;

class Function
{
public:
    double operator()(double _x);
    //double operator()(double _x, int n);
    double diff(double _x);
};

struct Point
{
    Point(double x, double y, double z)
    {
        this->x = x;
        this->y = y;
        this->z = z;
    }
    Point(){}
    double x;
    double y;
    double z;
    int k = -1;
};

vector<Point> *vec = NULL;

void init(Point points[], int k)
{
    vec = new vector<Point>();
    for(int i = 0; i < k; i++)
    {
        points[i].k = i;
        vec->push_back(points[i]);
    }
}


// double gdq(int k)//差商值，不带导数
// {
//     vector<Point>::iterator it;
//     if(k == 0)
//     {
//         it = vec->begin();
//         return (*it).y;
//     }
//     else
//     {
//         double result = 0;
//         for(it = vec->begin(); it != vec->end(); it++)
//         {
//             int j = (*it).k;
//             if(j > k)
//             {
//                 return result;
//                 break;
//             }
//             double a = (*it).y;
//             double b = 1.0;
//             for(int i = 0; i <= k; i++)
//             {
//                 if(i != j)
//                 {
//                     b = b*((*it).x-(*vec)[i].x);
//                 }
//             }
//             result = result + a/b;
//         }
//         return result;
//     }
// }

double h_gdq(int t, int n)//k从1开始,n要多加1
{
    double p[n][n+2];
    vector<Point>::iterator it;

        for(int i = 0; i < n; i++)
        {
            for(it = vec->begin(); it != vec->end(); it++)
            {
                if(i == (*it).k)
                {
                    p[i][2] = (*it).y;
                    p[i][1] = (*it).z;
                    p[i][0] = (*it).x;
                }
            }
        }//将第三列设置为函数值,第二列为导数值，第一列设置为x值
        for(int m = 0; m < t-1; m++)
        {
            for(int i = m; i < t; i++)
            {
                if(p[i+1][0] == p[i][0] && m == 0)
                {
                    p[i+1][3] = p[i+1][1];//若前后两个相等，设置为导数值
                }
                else
                {
                    p[i+1][m+3] = 1.0*(p[i+1][m+2]-p[i][m+2])/(p[i+1][0]-p[i-m][0]);
                }
            }
        }
    return p[t-1][t+1];
}



double gcm(double x,double k)//保存连乘结果
{
    double result = 1;
    vector<Point>::iterator it;
    int count = 0;
    for(it = vec->begin(); it != vec->end(); it++)
    {
        if(count ==k )
        {
            break;
        }
        result = result * (x-(*it).x);
        count++;
    }
    return result;
}

double Newton_y(double x,int n)//计算插值结果
{
    vector<Point>::iterator it;
    double result = 0;
    int index = 0;
    for(it = vec->begin(); it != vec->end(); it++)
    {
        if(index == 0)
        {
            result = result+((*it).y);
        }
        else
        {
            int k = (*it).k;
            double a = h_gdq(k+1,n+2);
            double b = gcm(x,k);
            result = result+a*b;
        }
        index++;
    }
    return result;
}